# add all imports, including for SQ2 and 3
import pandas as pd
import numpy as np

import plotly.io as pio
import plotly.express as px
import plotly.graph_objects as go
from plotly.offline import init_notebook_mode

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller

import warnings

init_notebook_mode(connected=True)
pio.renderers.default = "plotly_mimetype+notebook"
warnings.filterwarnings('ignore')

# Load data
def load_data(file_path):
    """
    Loads the data from the provided Excel file path, skipping rows that do not contain actual data.

    Parameters:
        file_path (str): Path to the Excel file containing the data.

    Returns:
        A dataframe containing the data from the Excel file.
    """
    df = pd.read_excel(file_path, sheet_name='Sheet 1', skiprows=10, index_col=0)
    df.dropna(how='any', inplace=True)
    return df

# Load Dataset
file_path = r"data/SQ1/Freight&Mail_all_countries_2003_quarterly.xlsx"
df_air = load_data(file_path)
# Create a Plotly figure
fig = go.Figure()
x = df_air.columns  
y_labels = df_air.index  
# Loop through each country and add a trace for it, excluding the United Kingdom
for country in y_labels:
    if country == "United Kingdom":
        continue
    # Replace ":" with 0 for missing data
    y_values = df_air.loc[country].replace(":", 0).values
    fig.add_trace(go.Scatter(x=x, y=y_values, mode='lines', name=country))

# Update layout with titles and labels
fig.update_layout(
    title="Freight and Mail Air Transport by Country Over the Years",
    xaxis_title="Year",
    yaxis_title="Freight and Mail (tonnes)",
    legend=dict(
        title="Country",
        x=1.05,
        y=1,
        font=dict(size=10)
    )
)

# Rotate x-axis labels and adjust layout for better readability
fig.update_xaxes(tickangle=45)
fig.update_layout(margin=dict(t=40, r=150, b=40, l=40))

# Show the plot
fig.show()

# Load Dataset
file_path = r"data/SQ1/Freight&Mail_five_countries_2003_quarterly.xlsx"
df_air = load_data(file_path)

# Preprocess data
def preprocess_data(df, country = 'Germany'):
    """
    Divides the dataframe into two lists. Train holds data from before 2020 Q1, 
    and test holds data from and after 2020 Q1.

    Parameters:
        df (pandas.DataFrame): dataframe of the data.
        country (str, optional): The country to process. Default is "Germany".

    Returns:
        A pandas.Series with train data and a list with test data.
    """
    time_series = df.T.loc[:, country]  # Germany for default
    train = time_series[(time_series.index < '2020-Q1')]  # training set
    # val = time_series[(time_series.index < '2018-Q1')]
    test = time_series[time_series.index >= '2020-Q1']  # test
    return train, test

def check_stationarity(train, threshold = 0.05):
    """
    Checks if a time series is stationary by applying the Augmented Dickey-Fuller (ADF) test and, 
    if necessary, the series is differenced to achieve stationarity. This step is critical for ARIMA 
    models, which require stationary data to provide reliable forecasts.

    Parameters:
        train (pandas.Series): The time series data to be tested and, if needed, differenced.
        threshold (float, optional): The p-value threshold for stationarity. If the ADF test p-value is above this level, 
            differencing will be applied to make the series stationary. Default is 0.05.
            
    Returns:
        The number of differences (d) needed to achieve stationarity, where 0 indicates 
        that the series is already stationary.    
    """
    # Perform ADF test on the original series
    result = adfuller(train.dropna())  # Ensure no NaN values
    print("ADF Test on original training data")
    print("p-value:", round(result[1], 4))
    d = 0
    # Check if differencing is needed (if p-value > 0.05, indicating non-stationarity)
    while result[1] > threshold:
        # Perform ADF test on the differenced series
        train = train.diff().dropna()
        result = adfuller(train)
        print("\nADF Test on differenced training data")
        print("p-value:", round(result[1], 4))
        d += 1
    else:
        print(f"The null hypothesis is rejected, which indicates that the training data is stationary and does not require differencing. Return parameter d = {d}")
        return d

def find_best_combination(train, d, display_models = True):
    """
    Identifies the optimal ARIMA or seasonal ARIMA model parameters for a given time series dataset 
    by iterating over combinations of non-seasonal (p, q) and seasonal (P, Q) parameters. The 
    best model is selected based on the Akaike Information Criterion (AIC). Lower AIC is better.

    Parameters:
        train (pandas.Series): The time series data for training the model.
        d (int) The degree of differencing required for stationarity, previously determined.
        display_models (bool, optional): If True, displays a summary of the top models with the lowest AIC scores. Default is True.

    Returns:
        The ARIMA or seasonal ARIMA model with the lowest AIC value.

    Notes:
        Seasonal period is fixed to 4 (for quarterly data).
    """
    # Convert the training data to an array
    train_array = np.asarray(train, dtype=float)

    # Define ranges for non-seasonal p and q
    non_seasonal_p_values = range(1, 4)  # Adjust as needed
    non_seasonal_q_values = range(1, 4)  # Adjust as needed

    # Define ranges for seasonal p and q
    seasonal_p_values = range(1, 4)  # Adjust as needed
    seasonal_q_values = range(1, 4)  # Adjust as needed

    # Initialize lists to store results
    non_seasonal_results = []
    seasonal_results = []
    # Loop over combinations of non-seasonal p and q values
    print("\nStart looping over combinations\n")
    for p in non_seasonal_p_values:
        for q in non_seasonal_q_values:
            model_nonseasonal = ARIMA(train_array, order=(p, d, q))
            result_nonseasonal = model_nonseasonal.fit()
            non_seasonal_results.append((p, q, result_nonseasonal.aic, result_nonseasonal))
            # Loop over combinations of seasonal p and q values
            for sp in seasonal_p_values:
                for sq in seasonal_q_values:
                    model_seasonal = ARIMA(train_array, order=(p, d, q), seasonal_order=(sp, d, sq, 4))
                    result_seasonal = model_seasonal.fit()
                    seasonal_results.append((p, q, result_seasonal.aic, sp, sq, result_seasonal))
    print("Loop over combinations done\n")
    # Combine results and find the model with the smallest AIC
    all_results = non_seasonal_results + seasonal_results
    best_model_info = min(all_results, key=lambda x: x[2])  # Get the model with the smallest AIC
    best_model_aic = best_model_info[2]
    best_model = best_model_info[-1]
    # Determine whether the best model is seasonal or non-seasonal
    if len(best_model_info) == 4:
        model_type = "Non-seasonal model"
    else:
        model_type = "Seasonal model"
    print(f"The best model is {model_type} p={best_model_info[0]}, q={best_model_info[1]}, P={best_model_info[3]}, Q={best_model_info[4]} and AIC={best_model_aic}\n")
    if display_models == True:
        # Display the training results
        non_seasonal_df = pd.DataFrame(non_seasonal_results, columns=['p', 'q', 'AIC', 'Model']).sort_values(by='AIC')
        seasonal_df = pd.DataFrame(seasonal_results, columns=['p', 'q', 'AIC', 'P', 'Q', 'Model']).sort_values(by='AIC')
        # Display the training parameters
        print("Non-seasonal models:")
        print(non_seasonal_df.loc[:, ['p', 'q', 'AIC']].head(3))

        print("\nSeasonal model:")
        print(seasonal_df.loc[:, ['p', 'q', 'AIC', 'P', 'Q']].head(3))
    else:
        pass
    
    return best_model

def forecast_plot(best_model, train, test):
    """
    Generates a forecast for the test period using the best-fitting seasonal ARIMA model and plots 
    the training data, actual test data, and forecasted values. This visualization aims to model a 
    "no-pandemic" scenario by forecasting based on pre-pandemic data and assessing recovery time by 
    identifying when actual volumes return close to expected levels.

    Parameters:
        best_model (statsmodels.tsa.arima.model.ARIMAResults): The optimized ARIMA model fitted on the training data.
        train (pandas.Series): Historical time series data used to train the model, ending in 2019.
        test (pandas.Series): Time series data for the test period (2020-2024) to compare against the forecast.

    Returns:
        None
        Displays an interactive plot of the training data, forecast, and test data. Also, prints the 
        recovery time, defined as the time point where the actual data is closest to the forecasted 
        values, suggesting alignment with pre-pandemic trends.

    Notes:
        The function computes the absolute difference between the forecast and actual values for each 
        point in the test period to determine the "recovery time"—the point when the volume aligns most 
        closely with the forecast, indicating a return to expected transport levels.
    """
    train_array = np.asarray(train, dtype=float)
    test_array = np.asarray(test, dtype=float)
    
    # Forecast future values using the best model
    test_array = np.insert(test_array, 0, train_array[-1])
    forecast_index = np.insert(test.index.values, 0, train.index[-1])  # Set forecast index to align with test data
    forecast_steps = len(test)  # Number of steps to forecast, matching the test set length
    forecast = best_model.get_forecast(steps=forecast_steps)
    forecast_mean = forecast.predicted_mean  # Forecasted mean values
    forecast_mean = np.insert(forecast_mean, 0, train_array[-1])
    # forecast_ci = pd.DataFrame(forecast.conf_int(), index=forecast_index[0:-1], columns=['lower', 'upper'])  # Confidence intervals
    
    # Calculate yearly difference and cumulative sum
    absolute_difference  = np.abs(test_array[2:] - forecast_mean[2:])  # Calculate yearly difference excluding first aligned value

    # Find the index of the smallest difference
    recovery_index = np.argmin(absolute_difference)

    # Get the corresponding date for the closest recovery time
    recovery_time = forecast_index[recovery_index + 2]  # +1 to account for the shifted index

    # Display the closest recovery time and its corresponding difference
    print("\nClosest Recovery Time:", recovery_time)
    print(f"Delay at Closest Recovery Time:{round(absolute_difference[recovery_index], 3)} (tonnes)\n")
    # Plot     
    # Create the figure
    fig = go.Figure()

    # Add the training data line
    fig.add_trace(go.Scatter(
        x=train.index,
        y=train_array,
        mode='lines',
        name='Training Data',
        line=dict(color='royalblue')
    ))

    # Add the actual data line
    fig.add_trace(go.Scatter(
        x=forecast_index,
        y=test_array,
        mode='lines',
        name='Actual Data',
        line=dict(color='orangered')
    ))

    # Add the forecast line
    fig.add_trace(go.Scatter(
        x=forecast_index,
        y=forecast_mean,
        mode='lines',
        name='Forecast',
        line=dict(color='seagreen')
    ))

    # Add the confidence interval as a shaded area
    # fig.add_trace(go.Scatter(
    #     x=forecast_index[0:-1].tolist() + forecast_index[0:-1][::-1].tolist(),
    #     y=forecast_ci['upper'].tolist() + forecast_ci['lower'][::-1].tolist(),
    #     fill='toself',
    #     fillcolor='rgba(0, 0, 255, 0.2)',
    #     line=dict(color='rgba(255, 255, 255, 0)'),
    #     hoverinfo="skip",
    #     showlegend=True,
    #     name='Confidence Interval'
    # ))

    # Update layout
    fig.update_layout(
        title=f'Best ARIMA Model Forecast of Air Transport Volume in {train.name}',
        xaxis_title='Time',
        yaxis_title='Volume of Goods Transported (tonnes)',
        legend=dict(x=0.01, y=0.99),
        xaxis=dict(tickangle=45)
    )

    # Show the figure
    fig.show()

train, test = preprocess_data(df_air, country = "Germany")
d = check_stationarity(train)
best_model = find_best_combination(train, d, display_models = True)
forecast_plot(best_model, train, test)

ADF Test on original training data
p-value: 0.3957

ADF Test on differenced training data
p-value: 0.0023
The null hypothesis is rejected, which indicates that the training data is stationary and does not require differencing. Return parameter d = 1

Start looping over combinations

Loop over combinations done

The best model is Seasonal model p=2, q=2, P=1, Q=1 and AIC=1518.8166590460717

Non-seasonal models:
   p  q          AIC
0  1  1  1655.715411
3  2  1  1659.054358
6  3  1  1659.671857

Seasonal model:
    p  q          AIC  P  Q
36  2  2  1518.816659  1  1
0   1  1  1519.620156  1  1
4   1  1  1519.807188  2  2

Closest Recovery Time: 2020-Q3
Delay at Closest Recovery Time:8927.282 (tonnes)

train, test = preprocess_data(df_air, country = "Netherlands")
d = check_stationarity(train)
best_model = find_best_combination(train, d, display_models = False)
forecast_plot(best_model, train, test)

ADF Test on original training data
p-value: 0.5825

ADF Test on differenced training data
p-value: 0.0027
The null hypothesis is rejected, which indicates that the training data is stationary and does not require differencing. Return parameter d = 1

Start looping over combinations

Loop over combinations done

The best model is Seasonal model p=1, q=1, P=2, Q=2 and AIC=1431.978071528646


Closest Recovery Time: 2020-Q3
Delay at Closest Recovery Time:4105.561 (tonnes)

train, test = preprocess_data(df_air, country = "France")
d = check_stationarity(train)
best_model = find_best_combination(train, d, display_models = False)
forecast_plot(best_model, train, test)

ADF Test on original training data
p-value: 0.706

ADF Test on differenced training data
p-value: 0.0038
The null hypothesis is rejected, which indicates that the training data is stationary and does not require differencing. Return parameter d = 1

Start looping over combinations

Loop over combinations done

The best model is Seasonal model p=1, q=1, P=1, Q=1 and AIC=1508.8928874466699


Closest Recovery Time: 2022-Q3
Delay at Closest Recovery Time:567.963 (tonnes)

train, test = preprocess_data(df_air, country = "Belgium")
d = check_stationarity(train)
best_model = find_best_combination(train, d, display_models = False)
forecast_plot(best_model, train, test)

ADF Test on original training data
p-value: 0.6773

ADF Test on differenced training data
p-value: 0.0186
The null hypothesis is rejected, which indicates that the training data is stationary and does not require differencing. Return parameter d = 1

Start looping over combinations

Loop over combinations done

The best model is Seasonal model p=2, q=3, P=2, Q=1 and AIC=1459.4987097091096


Closest Recovery Time: 2022-Q1
Delay at Closest Recovery Time:15160.374 (tonnes)

train, test = preprocess_data(df_air, country = "Italy")
d = check_stationarity(train)
best_model = find_best_combination(train, d, display_models = False)
forecast_plot(best_model, train, test)

ADF Test on original training data
p-value: 0.959

ADF Test on differenced training data
p-value: 0.2826

ADF Test on differenced training data
p-value: 0.0
The null hypothesis is rejected, which indicates that the training data is stationary and does not require differencing. Return parameter d = 2

Start looping over combinations

Loop over combinations done

The best model is Seasonal model p=2, q=1, P=1, Q=1 and AIC=1272.1595460783876


Closest Recovery Time: 2021-Q1
Delay at Closest Recovery Time:73414.697 (tonnes)

import pandas as pd
import plotly.graph_objects as go

def load_and_display_data(df_path, annotation_text, title_template_total, y_axis_title):
    """
    Loads the data from the provided Excel file path, displays the first few rows,
    and visualizes the data for import and export volumes over time per country using Plotly.
    Additionally, calculates and plots the total import and export values across all countries.

    Parameters:
        df_path (str): Path to the Excel file containing the trade data.
        title_template_total (str): Template for the plot title, including placeholders for {country}.
        annotation_text (str): Text for annotation (e.g., "COVID-19").
        y_axis_title (str): Label for the y-axis.

    Returns:
        None
    """
    # Load the data from the Excel file
    trade_data = pd.read_excel(df_path)

    # Pivot the data to organize it for plotting
    pivot_data = trade_data.pivot_table(values='primaryValue',
                                        index=['period', 'reporterDesc'],
                                        columns='flowDesc',
                                        aggfunc='sum').reset_index()


    # Calculating total import and export values across all countries for each year
    total_trade = pivot_data.groupby('period')[['Import', 'Export']].sum().reset_index()

    fig = go.Figure()
    fig.add_trace(go.Scatter(x=total_trade['period'], y=total_trade['Import'],
                             mode='lines+markers', name='Total Import'))
    fig.add_trace(go.Scatter(x=total_trade['period'], y=total_trade['Export'],
                             mode='lines+markers', name='Total Export'))

    fig.update_layout(
        title=title_template_total,
        xaxis_title="Year",
        yaxis_title=y_axis_title,
        legend_title="Trade Flow",
        xaxis=dict(tickmode='linear'),
        shapes=[dict(
            type='line', x0=2019, x1=2019, y0=0,
            y1=max(total_trade['Import'].max(), total_trade['Export'].max()),
            line=dict(color='red', dash='dash')
        )]
    )
    fig.add_annotation(x=2019, y=max(total_trade['Import'].max(), total_trade['Export'].max()),
                       text=annotation_text, showarrow=True, arrowhead=1, ax=-40, ay=-30)

    fig.show()

df_path = r'data/SQ2/Sorted_TradeData-27.xlsx'
load_and_display_data(df_path, "COVID-19", "Total Oil and Gas Trade Data Across All Countries (Import and Export)" ,"Total Trade Value ($ Dollar)")

df_path = r'data/SQ2/Sorted_TradeData_30.xlsx'
load_and_display_data(df_path, "COVID-19", "Total Pharmaceutical Products Trade Data Across All Countries (Import and Export)", "Trade Value ($ Dollar)")

df_path = r'data/SQ2/Sorted_TradeData_85.xlsx'
load_and_display_data(df_path, "COVID-19", "Total Electrical Machinery and Equipment Trade Data Across All Countries (Import and Export)", "Trade Value ($ Dollar)")

df_path = r'data/SQ2/TradeData020708.xlsx'
load_and_display_data(df_path, "COVID-19", "Total Perishables Trade Data Across All Countries (Import and Export)", "Total Trade Value ($ Dollar)")

file_path = "data/SQ3/Germany_Air_Trade_Data.xlsx"
file_path_sea = "data/SQ3/Germany_Sea_Trade_Data.xlsx"
df_germany = pd.read_excel(file_path)
df_germany_sea = pd.read_excel(file_path_sea)


df_germany_m = df_germany[df_germany["Trade Flow"] == "M"].copy()
df_germany_sea_m = df_germany_sea[df_germany_sea["Trade Flow"] == "M"].copy()


df_germany_m["Trade Value (US$)"] = df_germany_m["Trade Value (US$)"].replace('[\\$,]', '', regex=True).astype(float)
df_germany_sea_m["Trade Value (US$)"] = df_germany_sea_m["Trade Value (US$)"].replace('[\\$,]', '', regex=True).astype(float)


df_germany_m["Transport Mode"] = "Air"
df_germany_sea_m["Transport Mode"] = "Sea"
df_combined = pd.concat([df_germany_m, df_germany_sea_m])


fig = px.line(
    df_combined, 
    x="Period", 
    y="Trade Value (US$)", 
    color="Transport Mode", 
    title="Trade Flow (Import) for Germany by Year",
    labels={"Period": "Year", "Trade Value (US$)": "Trade Value (US$)"}
)

fig.show()

file_path = "data/SQ3/Germany_Air_Trade_Data.xlsx"
file_path_sea = "data/SQ3/Germany_Sea_Trade_Data.xlsx"
df_germany = pd.read_excel(file_path)
df_germany_sea = pd.read_excel(file_path_sea)


df_germany_x = df_germany[df_germany["Trade Flow"] == "X"].copy()
df_germany_sea_x = df_germany_sea[df_germany_sea["Trade Flow"] == "X"].copy()


df_germany_x["Trade Value (US$)"] = df_germany_x["Trade Value (US$)"].replace('[\\$,]', '', regex=True).astype(float)
df_germany_sea_x["Trade Value (US$)"] = df_germany_sea_x["Trade Value (US$)"].replace('[\\$,]', '', regex=True).astype(float)


df_germany_x["Transport Mode"] = "Air"
df_germany_sea_x["Transport Mode"] = "Sea"
df_combined = pd.concat([df_germany_x, df_germany_sea_x])


fig = px.line(
    df_combined, 
    x="Period", 
    y="Trade Value (US$)", 
    color="Transport Mode", 
    title="Trade Flow (Export) for Germany by Year",
    labels={"Period": "Year", "Trade Value (US$)": "Trade Value (US$)"}
)

fig.show()

Project Group - 8¶

Members & Student numbers:¶

Introduction¶

Research Topic:¶

Research Objective¶

SMART Criteria¶

Main Research Question:¶

Data¶

Subquestion 1¶

Data Pipeline 1¶

Data Visualisation¶

Step 1. Preprocess data¶

Step 2. Stationary Check¶

Step 3. Building the seasonal ARIMA Model¶

Step 4. Forecast and Plot¶

Step 5. Recovery time¶

Germany¶

Netherlands¶

France¶

Belgium¶

Italy¶

Subquestion 2¶

Data Pipeline 2¶

Subquestion 3¶

Data Pipeline 3¶

Conclusion¶

Discussion¶

Contribution Statement¶